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Noise and bias in square-root compression schemes 
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ABSTRACT 

We investigate data compression schemes for proposed aU-sky diffraction-hmited 
visible/NIR sky surveys aimed at the dark energy problem. We show that lossy square- 
root compression to 1 bit of noise per pixel, followed by standard lossless compression 
algorithms, reduces the images to 2.5-4 bits per pixel, depending primarily upon the 
level of cosmic-ray contamination of the images. Compression to this level adds noise 
equivalent to < 10% penalty in observing time. We derive an analytic correction to 
flux biases inherent to the square-root compression scheme. Numerical tests on simple 
galaxy models confirm that galaxy fluxes and shapes are measured with systematic 
biases < 10~^ induced by the compression scheme, well below the requirements of 
supernova and weak gravitational lensing dark-energy experiments. An accompanying 
paper ( Vanderveld et a/.|2009 ) bounds the shape biases using realistic simulated images 



of the high-Galactic-latitude sky. The square-root preprocessing step has advantages 
over simple (linear) decimation when there are many bright objects or cosmic rays in 
the field, or when the background level will vary. 

Subject headings: Data Analysis and Techniques 



1. Introduction 

The unexpected discovery of the acceleration of the Hubble expansion of the Universe has 
inspired numerous proposals for large-scale observational projects to seek further constraints on 
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this "dark energy" problem (see Frieman et al. 2008, for a review of dark energy). Several such 
proposals involve space-based visible or NIR imaging surveys of a substantial fraction of the full 
celestial sphere at resolution near the diffraction limit of 1-2 meter telescopes. Such surveys would 
also be extraordinarily valuable for general astrophysics, but the amount of data involved is large. 
The number of pixels in a survey of a fraction /gky of the sky, with pixel scale p and A^exp exposures 
per sky location, will be 

,14 ( J^Y"^ /skyiVexp 
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pix 



2.7 X 10' 
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A crate filled with disk drives is all that is needed these days for transmission of hundreds of 
terabytes of data between terrestrial sites. The transmission of this data from a spacecraft at the 
Earth-Sun L2 Lagrange point 1.5 million km away is, however, a substantial expense. If the data 
require an average of bpp bits per pixel to transmit, with a telemetry rate r, and a survey duration 
of T, then the number of hours of downlink per day required will be 

-1 



downlink time = bpp x 0.33 
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day 1.25/yr V^^OMbps^ 

To avoid excessive resource use both on the spacecraft and ground stations, the data volume must 
be reduced well below the 16 bits per pixel that is typical for digitization of astronomical array 
detectors. In this paper we investigate means for compressing astronomical imaging data well below 
this, ideally bpp < 4. 

Significant compression of the detector output signals will require a lossy compression step 
to reduce the bit depth of the noise — sometimes called "pre-processing" — followed by a lossless 
algorithm to reduce the mean bit depth to the signal entropy level. Because most of the pixels in 
an extragalactic visible/NIR exposure will be close to the sky level, it is clear that the primary 
determinant of the final bpp value will be the number of noise bits retained by the lossy step in pixels 
near sky level — extensive demonstration of this for ground-based images is done by [Pence et al. 



( 2009 ) . Optical astronomers are averse to lossy compression for fear of losing the information from 



their dearly acquired photons or fear of inducing biases in high-precision measurements. It should 
be pointed out, however, that lossy compression is already being performed on essentially all data 
when the analog detector outputs are digitized. Analog-to-digital gains are usually set, however, 
so that the errors induced by digitization are below the read noise of the detector system. In radio 
astronomy it is well known that digitization to only 1 or 2 bits per sample can capture nearly all the 
signal information when S/N per sample is low ( Weinreb|1963 Thompson et a/.|2001 ). The critical 
parameter for information loss is the number of bits of noise retained, i.e. the ratio b = a]\f/A 
of the signal noise to the digitization step size. We seek a lossy compression scheme which holds 
b fixed over the full dynamic range of the detector. When the signal noise is dominated by shot 
noise from the detected photoelectrons, it is straightforward to show that a square-root compression 
scheme holds b constant. The application of square-root compression to Poisson-dominated data 



has been discussed by Gowen & Smith (2003) and is part of the signal chain for the James Webb 



Space Telescope ( Nieto-Santisteban et al.l|1999 ) and Supernova Acceleration Probe (Lin Sz Marriner 
2003 ) spacecraft designs. Further discussion is in [Seaman, White, k Pence] ( [20091 ) . In ^ we review 



- 3 - 



the algorithms of such a scheme and ^ shows that the square-root compression algorithm increases 
the noise by a fraction 1/246^ for cases of interest. 

Once a lossy compression algorithm is determined to reduce the data volume with acceptably 
small degradation of the noise level, it is critical to determine whether the codec process will induce 
any bias on the astronomical measurements from the data. Future dark energy experiments will 
require very demanding, high-precision measurements. For use of Type la supernovae as distance 
indicators, systematic errors in flux measurements must be <C 1%. For determination of galaxy 
redshifts using the "photo-z" method, similarly small systematic errors must be present in galaxy 
photometry. We target codec-induced systematic errors in flux to be < 10~^ so that they contribute 
negligibly to the photometric error budget. 

Even more demanding are programs for weak gravitational lensing (WL) measurements, which 
depend upon determination of the shapes (ellipticities) of > 10^ galaxies on the sky. Huterer et 



al. (2006) and Amara &: Refregier ( 2008| ) both conclude that measurements of galaxy ellipticities 
must have systematic errors of < 10~^ in order to remain subdominant to statistical errors in the 
WL experiments. Because there are many possible sources of shape measurement error — primarily 
uncorrected instrumental image distortions — we would like to have shape biases induced by the 
codec be a minor contributor to the shape error budget, i.e. the shift Jcj in each component Cj of 
galaxy ellipticity should on average be < 10~'*ej. The goal of this paper is to examine the biases 
induced by square-root compression down to 5 = 1 bit of noise per pixel, and to show how these 



biases can be reduced to negligible levels. A companion paper (Vanderveld et al. 2009) tests for 



codec-induced shape bias on realistic simulations of space-based CCD imaging. 

Many lossy image compression algorithms are based on decomposition of 2d images over a 
function library, followed by truncation of the terms in the library — for example wavelet transforms 
( White et al.|[l992||Press|1992p or the JPEG algorithms. Even more radical compression is possible 



with "compressed sensing" techniques if the data are known to be sparse in some library (Candes, 



Romberg, &; T. Tao 2006). For WL applications we are wary of lossy compression algorithms 
that examine groups of pixels and might hence induce correlated codec errors over multiple pixels. 
We therefore focus here on a lossy step that compresses data strictly on a pixel- by-pixel basis. If 
single-pixel codec errors are unbiased and completely uncorrelated between pixels, object shape 
measurements should be unbiased. We investigate the single-pixel codec in this work, leaving open 
the question of whether grouped-pixel algorithms can meet the shape bias requirements of weak 
lensing surveys while obtaining better compression than single-pixel codecs. 

In ^we show that there is a significant bias in the mean values reconstructed by the square-root 
class of codecs. We also show, however, that a straightforward adjustment of the decoded values 
can remove this bias to remarkable accuracy, largely independent of the details of the probability 
distribution of the raw data. We therefore expect the codec to produce unbiased flux and shape 
measurements, which we test with simplified galaxy images in ^ [Vanderveld et al.\ ( |2009D test for 
shape-measurement bias on data that is a realistic simulation of images expected from the Joint 
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standard lossless algorithms, and measure the number of bits per pixel required to transmit the 
fully compressed images. In the conclusion we summarize our findings, namely that it is possible 
to compress JDEM-like images to 3-4 bits per pixel with only 4% penalty in image noise, and to 
reconstruct them such that biases in flux and shape measurements induced by the codec are not 
detected, down to « 10"'* levels. 



We consider most generally a compression algorithm that maps the digitized detector output 
value into a compressed code i whenever Ni — Ai/2 < N < Ni + Aj/2. Here Ni and Aj are 
the center and width, respectively, of the data range that codes to compressed value i. Note that 
Ni-^-i — Ni = (Aj + Aj4-i)/2. We will assume that is a debiased value that is linearly related to 
the estimated number of photoelectrons e by = e/g, where g is the (inverse) gain of the system. 
We will assume initially that decompression consists of mapping i ^ Ni, i.e. we restore code i to 
the mean of the input values that it codes. In a later section we will suggest an alteration to the 
decompressed values which reduces the bias in reconstructed values. 

Any vector of increasing Ni (or equivalently a vector of Aj > 0) defines a codec algorithm. We 
expect the bias and noise characteristics of the codec to depend primarily on the the number of 
bits of noise retained in the compression, which we define as 



where cr at is the RMS noise level of a pixel and A is the typical Aj for the codes i that span the 
noise distribution of the pixel. Many astronomical detectors have noise described by 



with some read noise level RN added to the Poisson variance of the photoelectrons. In this case, 
if we wish to keep b independent of the signal level, we require 



2. Compression algorithm 



b = ctat/A, 



(3) 





(4) 



A 



di 



UN/b 



(5) 



dN 



di 
dN 



b 



(6) 



ViN + c)/g 

^/k{N + c), k = Agb^. 



(7) 



' http : // j dem . gsf c . nasa . gov 



- 5 - 



Since the compressed code i must be integral, we can more precisely express the nominal square-root 
compression algorithm as 

' ' (8) 



^/k{N + c) = 2bVe + RN^ 



where, in this equation only, [x] is the nearest integer to x. Gowen Sz Smith (2003) discuss more 
general offsets under the square root, but we stick to these values since they hold o-jy/ A fixed across 
the dynamic range. 



2.1. Range step sequences 

If we define 

= A, - Ai_i (9) 
then from Equation ([T]) we find that the square-root codec has 

2y^k{N + c) _ 2i 



k 



(10) 



i.e. the width Aj of each code range increases at a constant rate. If the raw data N are integral, 
then these Aj and A^ must also be integral. If for example we desire 6=1 bit of noise and have 
gain g = 0.5 e/ADU, then A^ = 1 produces the desired compression algorithm. When l/2b'^g is 
non-integral, however, we will more generally obtain (A^) = l/2b'^g. In fact we can define a class 
of compression algorithms as follows: 

1. Choose a sequence of integers {si, S2, . ■ ■ , Sn} that will repeat cyclically to give the incremental 
range steps {A[, A'2, . . .} . 



2. Define the parameter b via 



6^ = — -. (12) 



Choose a range Aq for the first compression code that extends from = to = {RN) / gb. 
[The codec algorithm can be extended to small negative values of N if desired using the fixed 
range width Aq for small negative i values.] 



Then the codec defined by this range-step sequence will have nearly constant number b of bits per 
pixel noise level under the Poisson+read noise data model. The more nearly uniform the values 
Si are, the better will be the codec performance, as shown below. For example we can produce a 
codec with 6=1 with g = Ihy using code steps A^ = {0, 1, 0, 1, 0, 1, . . .}. We will refer to this as 
a "5 = 1, {0, 1} codec." 
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2.2. Data model 

When considering the compression to act upon integer data, we denote the probabiHty of the 
raw pixel value being as P/v- For notational convenience we will assume that there exists a 
continuous and differentiable function P{x) over x > which matches P/v at x = iV. For example 
the Poisson formula for P/v is easily extended to non-integral values with a gamma function. 

In the CCD data the Poisson signal has a Gaussian read noise added to it. We will assume 
the read noise has zero mean. Before digitization we can hence consider the pixel photoelectron 
count e to be a real number, with probability distribution P{e). In the limit of high gain g we can 
consider the detector output to be a real number e or A^; at finite gain, Pn is again defined only at 
integral A^, but we can assume there exists some analytic extension to the full real domain P(N). 

Most of our results are general to any well-behaved P{N) or P(e), and we will note when the 
assumption of Poisson or Poisson+Gaussian noise models are important. 



3. Codec noise 

A raw pixel value will incur an error 5N = Ni — N upon being encoded and decoded through 
code i. If the input values are uniformly distributed across the range of code i, then the RMS error 
induced by the codec is fXcodec = \/ \fVl. If we make the further simplifying assumptions that the 
coding error is uncorrelated with the raw A^ and that A is nearly constant for all i accessed by 
P{N), then the variance of the output is 



'out 



1 



(^N + (^codec = (^N [^ + ] ■ (13) 



The integration time required for a shot-noise-limited astronomical exposure to reach fixed signal- 
to-noise ratio (S/N) will be inflated by 



a% ~ ^ ' 1252 • 



1 + (14) 



If we want the codec to cause < 10% penalty in exposure time, we hence need 5 > 0.91. In this 
paper we will focus our attention on this question: does a square-root algorithm with b > 1 cause 
any bias or degradation in astronomical measurements beyond the expected < 4% increase in noise 
level? 

A floor on the bandwidth required to transmit real images is given by the bpp found for 
pure-noise images. The minimum bpp required to for lossless compression of a sequence of codes i 
with probabilities pi is the Shannon entropy H = —'^pilog2Pi- For nearly-Gaussian probability 
distributions digitized to b bits per a, the Shannon entropy is H ^ log2 -v/2vre + log2 b{e.g. 



Romeo 



et al. 1999). For 6 = 1 a numerical evaluation of the Shannon entropy shows that > 2.1 bpp will 



be required. We construct a 16-bit integer FITS image consisting of random noise with e = 40 and 
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RN = 4, apply square-root compression with 6 = 1, then apply lossless bzip2 or CCSDS 121B 
( CCSDS|1997 ) compression. The resulting image requires 2.4-2.5 bits per pixel, within 0.3-0.4 bpp 
of the Shannon limit. We can therefore consider 2.4 bits per pixel to be a floor on the telemetry 
required when using a pixel-by-pixel lossy compression step, subject to < 10% integration-time 
penalty from codec errors. We note that the "raw" data digitized at g = 1 would have b = 7.5 and 
require 2.9 more bits per pixel. We verify that with square-root compression enabled, the bpp for 
pure-noise images is largely independent of e, RN, and g, as expected. 



4. Reconstruction bias 

The expectation value of the raw data N = (N) should be unchanged after compression and 
decompression. We will denote by 6N the bias in the expectation value of the data after the codec 
operates. In practice this bias is not a problem if it is independent of N, since it is easily corrected. 
If 6N varies with N and/or with the details of the P{N) distribution, then the response of the 
detector effectively becomes nonlinear and/or ambiguous to correct, which can bias photometry 
and shape measurements. 

Any codec bias could of course be corrected if the raw distribution P/v is known at each pixel. 
But the mean level N will vary from pixel to pixel due to background gradients and sources, so 
there is no single Pn that applies to the entire image. We would prefer to debias the data by 
decompressing code i to Ni — 5i, where 6i is some small correction that is independent of the raw 
probability distribution. There is no guarantee that such a solution exists but we will show that a 
good approximation does exist for cases when P{N) spans multiple codes. 



4.1. Real- valued data 

We begin with the case in which the raw value N can be considered a continuous real variable. 
This case is instructive if not necessarily realistic. The center and width Ni and Aj of the range 
for each code i can be set to any positive real numbers. The codec bias can be expressed as 

5N = -} duuP{Ni + u), (15) 

where u = N — Ni is the distance from the center of code step i. If we can find a set of 6i such that 

/■A,/2 

y / du5iP{Ni + u) = 5N (16) 
i 

holds for any distribution P{N), then we will have a complete de-biasing prescription for the codec. 
When the distribution P{N) is wider than a code step width A, it can be approximated by 
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the Taylor expansion of P{N) about Ni: 

P{Ni + u) = Pi + P[u + P'lv?l2 + i^"'nV6 + 0(n/aJv)^ 



(17) 



with P[ = P'{Ni), etc. We will truncate the series at cubic order henceforth. The contribution of 
the neglected terms will scale as (A/ctat)^ = b~^. We will later perform a numerical check of the 



significance of the neglected terms. Substituting the Taylor expansion into Equation (15) gives 

5iV « - ^ A, [P/Af /12 + PrAf/m] . (18) 



If this sum over i were an integral, we would evaluate this quantity through integration by 
parts. We can perform an analogous discrete process by defining ^4-1-1/2 = P{^i^ Ai/2) and noting 
that to our chosen order of Taylor expansion, we have 



^iP'i — Pi+1/2 — Pi-1/2 



A3 

izi p'" 
24 ^ ■ 



Substituting this into Equation (18) yields 

"A? 



5N 



E 



A5 



i-1/2) 



720 



-P" 



E 



^A^P,_l/2 + A^+lP.+l/2) + ^Pl" 



12 



(19) 

(20) 
(21) 



Recalling that A^ = A, — Aj_i, it is clear first that the codec bias is very small if the step width A 
is constant (A' = 0). This would be the case if the codec were a simple decimation over the range 
of the data. This is the case for normal digitization of the data. 

The perfect square-root compressor of real-valued raw data Equation ([s]) has the unique prop- 
erty that A[ = A' is constant. In this special case we can simplify further and used the Taylor 
expansion one more time to get -Pj+1/2 + Pi-1/2 = 2Pj + AfPl'/4, and now we have 



5N = Y,^i 



A2 

Pi + -^p:' ] + 



720' 



-P" 



(22) 



Second, it is clear that the square-root codec bias will be nearly independent of P{N) because 
AjPj — )• 1 in the limit of small Aj, so we expect there to be bias 



5N 



A^ 



1 



(23) 



A constant bias is relatively benign, for example does not influence the appearance of a sky- 
subtracted image. It is also easily remedied by decompressing code i to value Ni — 5^^\ where 
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jC') = A'/6. After making this adjustment, the residual bias, to third order in the Taylor expansion 
of each code range, becomes 



6N = yAi 



72 ' 720 ' 



(24) 



Further application of the "integration by parts" technique shows that the leading terms in this 
residual bias scale as / dN A'P'{N), which is zero if A' is constant across A^. Hence we expect the 
reconstruction bias to be very close to the constant 1/126^5 for the reconstruction of real- valued 
raw data with 6 > 1, i.e. whenever the intrinsic noise spans multiple compression code steps. 

To summarize the key results: a simple decimation codec (Aj =constant) reproduces the mean 
of the input data with no bias, and the square-root codec (A^ = constant) biases the mean by 



the constant value in Equation ( 23 ) independent of the distribution P{N) of the raw data. These 
statements are very good approximations when b > 1 and the raw data are real-valued. 



4.2. Integer data 



Next we consider the case when the raw data N must be integral. In this case Aj and A[ 
must be integral, and Ni is integral or half-integral depending upon whether Aj is odd or even, 



respectively. The integrations over N in Equation (15) must be replaced by sums over the allowed 
integral values of N. We obtain 



where the sum runs over j G {-(Aj - l)/2,-(Ai - 3)/2,. . . , (A^ - 3)/2, (Aj - l)/2}. 

We smoothly extend the discrete probabilities to a continuous function P{N) that main- 
tains the normalization condition J dN P(N) = 1, and we again perform a 3rd-order Taylor expan- 
sion of P{N) about the center Ni of each code range. The bias can now be written 



5N 



E 



(26) 



The sums of powers of integers (or half-integers) can be found from well-known formulaej^ The 
result is 



6N 



E 



12 



f A? -I) P! 



3A,2 - 7 
120 



P' 



(27) 



For example, 



http : / /mathworld . wolfram . com/FaulhabersFormula . html 
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The discrete version of integration by parts transforms this into an expression identical to Equa- 



tion (21) up to small corrections in the P'" term: 



6N = Y^ 



12 



-1/2 



1/2 



(Af-l)(2Af + 7) 
120 



(28) 



For an integer-valued codec we cannot necessarily maintain constant A^ for an arbitrary 6; for 
example with g = 1 and 6 = 1, we must have (A') = 1/2. Hence we must choose adjustments 5i to 



each decompressed value that will best cancel the bias in Equation (28). After further manipulation 



it is possible to show that the leading-order terms in the bias are cancelled, independent of the 
form of Pat, if 



i+2 



A' 



A' + A'_ 



12 



A,_i 



48 

Ai+2 - 2Ai+i + 2Ai_i 



Ai_5 



12 



48 



(29) 
(30) 



To summarize: the compression algorithm is defined by a repeating integer sequence {s,} for the 



A' values, along with a gain g and read noise RN , as described in ^ 2.1 The compressed values will 
have h = {2g{si)y^l^ bits per a of input noise. On decompression, each compressed code i should 
be replaced with the value Ni — Si, with 6i given by Equation (pOl) . 



One should prefer the most uniform possible sequences {si}. A constant value (nominally 
A' = 1) will yield a "perfect" square-root compression with minimal reconstruction bias. An 
alternating sequence, i.e. {0, 1}, has a constant bias correction 5i, but slightly larger biases. Longer 
{si} sequences have more complicated and less precise bias corrections. 



4.3. Numerical validation 

Figure [T] shows the noise and bias introduced by several versions of the square-root codec. In 
each case the raw data are drawn from a Poisson distribution with the mean number of photoelec- 
trons plotted on the x axis, ranging from 40-lOOe. The data are taken to have Gaussian read noise 
with zero mean and 4e RMS. Each plotted point is the result of lO"^ trials from this distribution. 
Each trial produces an integral raw N value and a real- valued decompressed value Ni — Si. 

We examine 4 different codecs. The red, blue, and green curves are designed to produce 6 = 1 
bit of noise. The red curve is a "perfect" square-root codec, with A' = 1 and g = 0.5e per ADU. 
The blue and green curves have g = 1 and {si} sequences of (0, 1) and (0, 0, 1, 1), respectively. The 
magenta curve is a codec with g = 1 and step sequence (0, 0, 1), which retains b = 1.22 bits of noise. 

The bottom panel plots the mean gNi and range gAi of each code step i for each of the four 
codecs, expressed in units of photoelectrons. The top panel shows the RMS fluctuation of the 
decoded data, relative to the RMS fluctuation of the input data. All of the codecs inflate the RMS 
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Fig. 1. — Performance of four square-root codecs on data with Poisson noise plus 4e RMS Gaussian read 
noise. The bottom panel shows the center values Ni and the range widths A.; of each step for the four codes. 
The red (solid), green (dotted), and blue (long dash) codecs each produce 6 = 1 bit per noise crjv of the input 
value, with decreasing levels of smoothness in the derivative A' of the range width. The magenta (short 
dash) codec retains more noise, b — 1.22. The top panel shows the increase in data RMS due to the codec, 
which is independent of input level. The central panel shows the bias in the reconstructed data stream, with 
the noise level of this measure indicated by the black error bar. Bias levels are O(O.OOle), in some cases 
consistent with zero within measurement noise. See text for further detail. 
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by almost exactly the factor (1 + 1/126^)^/^ predicted in Equation (13). In particular the b = 1 
codecs inflate the noise by 4%. 

The central panel plots the mean difference between decompressed and raw data. The expected 
noise in this measurement from the 10'' trials is plotted in black at left. The "perfect" codec, in 
red, has zero bias within errors, \6N\ < 0.001. This bias is, remarkably, < 10~^ of the input values, 
< 10~^ of the typical noise and code range A =7-10, and also < 1% of the bias 1/126^ = 0.083 
that the codec would have generated without any correction. The {0, l}-sequence, in green, shows 
some significant residual bias at the ±0.002 level, and the {0, 0, 1, l}-sequence codec (blue) has bias 



of ~ ±0.004 electrons. While these biases are quite small, they do highlight that Equation (30) 
ignores some higher-order terms that make the bias dependent on P/v, and that it is preferable to 
use the smoothest possible step sequences. The {0, 0, 1} codec also has some bias at the 0.001 level, 
despite having slightly less compression and lower noise than the "perfect" 6 = 1 codec in red. 

These numerical tests insure that square-root compression at 6 = 1 will introduce no significant 
biases into photometric measurements, and increases in errors will be limited to 4%. 



5. Shape measurement biases 



Using simple simulated galaxies and a very basic shape-measurement routine, we demonstrate 
here that the codec does not induce any biases in the measurements of galaxy fluxes, sizes, or 
ellipticities at the part-in- 10^ level. In this work we analyze simple galaxies with specified size and 
signal-to-noise ratios which span a range that we might expect in real images. A test for codec 



with much more realistic images is presented by Vanderveld et al. ( 2009 ) . In that work the galaxies 



have sizes, appearances, fluxes, and shapes that accurately mimic those we might expect from a 
visible survey telescope in orbit. 



5.1. Procedure 

The simulated images are postage stamps constructed as follows: 

1. Start with an exponential-proflle galaxy, with intensity / oc e~^^'^^. The exponential disk is 
convolved with a Gaussian of width a = ro/10, to suppress high frequencies associated with 
the r = cusp of the exponential disk. 

2. This circular galaxy is then sheared to an ellipticity e = {a? — b'^)/{a? + 6^) = 0.1, 0.3, or 0.7 
(a and b are the major and minor axes of the galaxy). The major axis of the elliptical galaxy 
can have position angle (/> = 0, 22.5°, or 45° from the x axis. 

3. The galaxy is (de-)magnified until the resulting galaxy has half-light radius of ri/2 = 2, 4, or 
6 pixels. 
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4. The galaxy profile is shifted by a random fraction of a pixel in x and y, then convolved with 
the square pixel response function and sampled on a 1-pixel grid to yield a noiseless postage 
stamp image. 

5. The flux of the galaxy is adjusted until its total signal-to-noise ratio will be S/N = 15, 50, 
or 300 when noise is added in the following step. The S/N is defined as 

{S/Nf= ^ (31) 

iSpixels ' 

where li and ai are the signal and noise in each image pixel. 

6. A sky level of 50 electrons is added to each pixel's expected signal Ij. Then an integral 
"observed" photoelectron count for each pixel is determined by drawing from a Poisson dis- 
tribution with mean /j. Gaussian read noise with zero mean and a = 4 electrons is added to 
give a floating-point measurement value for each pixel. 

At this point two digitized versions of the postage-stamp image are produced. For the "codec" 
image, the photoelectron count /j is scaled by gain g = 0.5e/ADU, digitized, compressed with a 
sequence 1 perfect square-root codec, and then decompressed to a floating point image following 
the debiasing procedures above. Note that the codec retains 6=1 bits per noise sigma on each 
pixel. For the "raw" image, we add Gaussian noise with RMS of ai/\/T2b to each pixel, so that the 
noise level of the raw and codec images will be expected to be equal. Then we scale by the gain g 
and digitize without square-root preprocessing. 

Both the raw and codec versions of the postage stamp galaxy are measured using a basic 

as follows: the pixel 

data are flt to an elliptical Gaussian model, with the flux, centroid, size, and ellipticity of the 
Gaussian allowed to vary. The sky background is fixed at 50/ g = 100 ADU. There are hence 6 
degrees of freedom in galaxy model (and the Gaussian galaxy model is not an exact fit to the 
underlying exponential-profile galaxy). 

5.2. Results 

For each choice of ellipticity e, orientation 0, size rx/2 of galaxy, we create 1000 simulated 
galaxies at S/N = 300, 30,000 realizations with S/N = 50, and 300,000 reahzations with S/N = 15. 
To test for biases induced by the codec, we measure the mean difference ecodec — ^raw between the 
ellipticities measured on the codec and the raw images. Similarly we measure the mean difference 
in the (log of) flux and size measured for the galaxies. To see if the codec induces extra shape 
noise, we compare the RMS variation of the fitted e in the codec data set to the RMS variation in 
the raw fits. We also check for additional variance in fiux and size measurements. 



"adaptive moments" shape- measurement algorithm (Bernstein &: Jarvis||2002 
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Fig. 2. — The biases (left) and relative noise level (right) of measurements on codec images, relative 
to raw images, are plotted for several quantities measured on simulated galaxies. The details of the 
simulation and the measurement process are in the text. From top to bottom, we plot the bias in 
the (log of) galaxy flux, (log of) galaxy size, and then the two components of galaxy ellipticity. Red 
(outlined) points are for galaxies with S/N = 15, green (filled) have S/N = 50, and blue (shaded) 
have S/N = 300, with each point plotted at its input ellipticity. All data are fully consistent with 
the codec images yielding no bias or change in noise level from the raw images — at accuracies of 
±10"^ for biases and < 1% in noise levels. [Noise levels for S/N = 300 cases are not plotted because 
there are too few trials for an accurate measure of RMS noise. 1 
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Figure [2] presents the biases and noise levels on the ellipticity, size, and flux measurements of 
galaxies with = 2 pixels and </> = 0. The results for larger galaxies and for different ellipticity 
orientations are similar and not plotted here. We find no evidence for any codec-induced bias in 
size, flux, or ellipticity measurements at any S/N or ellipticity. The number of trials is sufficient 
to determine these biases to ~ ±10"^ for each configuration, approximately 10 times below the 



maximum permissible ellipticity measurement errors in a large WL survey (Amara & Refregier 



2008 Huterer et al. 2006). For these simplified galaxies and shape measurement techniques, we 
confirm that the square-root codec at 6 = 1 does not induce any significant bias in weak-lensing 
shape measurements. 

The S/N = 50 and S/N = 15 cases have enough trials to determine the measurement noise 
in e, size, and flux to <C 1%. We find the raw and codec images to have identical measurement 
noise levels to this accuracy, confirming that the noise added by the codec is accurately described 



by Equation (13). 



6. Compression factors 

6.1. Simulated scene 

We create simulated images from a space-based visible survey telescope to estimate the bits of 
telemetry required per pixel (bpp) needed to transmit the images after they have been compressed 
with the square-root algorithm followed by a lossless compression algorithm. The simulated images 
have three components: the expected cosmic scene; a uniform sky background; and bright linear 
features produced by cosmic rays in the detectors. Our test images are aimed at reproducing 
the appearance of images taken with an instrument resembling the CCD camera proposed for the 
JDEM/IDECS mission. The relevant characteristics assumed are: a 1.5-meter diameter telescope 
operating in an L2 Lissajous orbit, where the sky background is dominated by zodiacal light; a 
200-second integration through a bandpass in the red part of the visible spectrum, with bandwidth 



ln(Amax/Amin) = 0.3; a fully-depleted, 200- micron-thick p-channel CCD (Roe et al. 2007); a scale 
of 0'.'14 per pixel; and a point spread function (PSF) with EE50 radius 0'.'14, after including all 
optical, mechanical, and detector-induced sources of image blur. 

The cosmic scene is taken from deep moderate-latitude exposures taken with the ACS/WFC 
F606W filter aboard the Hubble Space Telescope in the course of a search for very faint outer-solar- 
system objects ( Bernstein et al.|[2004 ). These exposures combine many HST exposures to produce 



an image whose noise is far below that expected from a single JDEM/IDECS exposure as taken 
above. The scene is scaled so that a source of magnitude Iab = 31.09 produces 1 photoelectron. 
The HST image is blurred and repixelized to match the JDEM/IDECS specifications described 
above. 



Diffuse sky emission is added assuming high-ecliptic-latitude zodiacal background at L2. We 
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estimate ~ 45 electrons per pixel in the detector. 



6.2. Simulated cosmic rays 

Cosmic rays in space are isotropic and uniformly impinge on a sphere. From ESA's Space 
Environment Information System (SPENVIS|^ program, the cosmic ray rate is found to be approx- 
imately 5 particles cm^ over Att steradian. This is the average rate for Galactic cosmic rays, 
which are dominated by protons. Low energy particles (mainly solar protons) have little integrated 
effect on the total rate even without shielding. We ignore periods of enhanced low momentum solar 
protons during solar storms. These episodes are known in advance and data taking is avoided. 

We use SPENVIS to generate a momentum spectrum for a particular orbit (L2 in this case) at 
a particular solar epoch behind appropriate shielding (3 mm Al assumed here). The Monte Carlo 
simulation samples this spectrum and for each momentum, an average ionization rate (e//im) is 
calculated for the target material (Si in this case) according to a modified Bethe-Bloch formalism 



(Amsler et al. 2008, Section 27). For the data generated and analyzed here, a single particle 
momentum was used corresponding to minimizing protons — approximately 3 GeV/c corresponding 
to 70 e/fim in silicon. 

For a given exposure time, particle flux and detector area, the total number of cosmic rays is 
calculated. We generate this number of cosmic rays uniformly over 47r sr. Each ray is randomly 
placed over the surface of detector and tracked through the detector material. Tracing is done in 
uniform step lengths. For each step, a Landau distribution is sampled to determine the number 
of ions produced. Hard scatters and (5-rays are not simulated. For each ion generated in a step, a 
Gaussian lateral charge diffusion distribution is randomly sampled. The diffusion constant is based 
on a measured 4 /im RMS for a 200 /xm thick CCD operated at lOOV bias voltage. This is linearly 
scaled by ion conversion height above the pixel plane. The diffused charge is projected to the pixel 
plane where its position is discretized and added to any previous pixel charge. 

In the simulated 200 s exposure, 7% of the pixels receive > 5e of CR-generated charge. We 
also simulate an image with triple this dose of cosmic rays, which is simply the sum of 3 200s CR 
simulations. 



6.3. Compression results 

The scene, cosmic ray, and sky components are summed, and then Poisson noise and 4e RMS 
read noise are added to each pixel. The noisy images are then scaled by gain g = 0.5e/ADU and 
digitized, then compressed using a perfect square-root algorithm with sequence 1. The compressed 



' http : //www. spenvis . oma.be/spenvis/ 
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images are saved as 16-bit integer FITS-format images. A segment of a simulated image, before 
compression, is shown in Figure [3} 

The Shannon entropy of the square-root-compressed image, calculated by direct measure of 
the Pi of each code in the image, is 2.8 bits per pixel. Tripling the cosmic-ray dose on the exposure 
increases the entropy to 3.6 bits per pixel, suggesting that the normal dose of cosmic rays is 
responsible for 0.4 bits per pixel of entropy in the image. Thus most of the 0.7-bit increase in 
entropy over the pure-noise image (as measured in ^ is due to cosmic rays. 

There are many possible lossless compression algorithms to apply to the square-root-compressed 
FITS images. This choice will of course not affect the biases in the data, but will affect the com- 
pression ratio and the robustness of the telemetry stream. Applying the commonly available bzip2 
algorithm to the FITS files yields the following results: 

• The images simulating a 200-second exposure and 200 seconds' worth of cosmic-ray tracks 
occupy 3.1 bits per pixel after square-root and bzip2 compression. 

• The information content of the images is dominated by the cosmic-ray component. If we 
triple the number of cosmic rays (mimicking a 600-second exposure, even though the cosmic 
scene and sky background are still 200 seconds' worth), we get a pessimistic value of 4.0 bits 
per pixel in the compressed images. 



The bzip2 performance is again 0.3-0.4 bpp above the entropy of the image. 

Full- frame bzip2 compression is not very robust: a single-bit transmission error can lead to 



loss of a full image. A more robust algorithm is CCSDS 121B (CCSDS 1997), which has already 
been implemented on > 25 space missions. We find the file size after square-root pre-processing 
and CCSDS 121B lossless compression to be within 0.1 bits per pixel of the bzip2 results. 

A different lossless compression scheme, optimized toward efficient compression of the linear 
cosmic-ray tracks that dominate the 2d image, might yield better compression yet. If further 
compression were desired, the most effective approach would be to identify cosmic-ray-contaminated 
pixels on that satellite and replace them with a constant tag value before compression, so that their 
entropy is not transmitted to the ground. 



6.4. Comparison to linear preprocessing 

Implementation of square-root preprocessing is very simple, requiring just a fixed lookup table 
to convert ADC output values to compressed codes. An even simpler algorithm is linear prescaling, 
i.e. decimation, in which the code range A is fixed over the full dynamic range. As noted above, a 
decimation algorithm should be free of bias, just as the corrected square-root algorithm is measured 
to be free of bias. What are the relative merits of the two lossy codecs? 
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Fig. 3. — Section of an image used to test the square-root+lossless compression scheme. This 
simulates a 200-second CCD exposure with a 1.5-meter space telescope observing near the ecliptic 
pole from an L2 orbit. Most pixels are dominated by sky and read noise, and celestial objects are 
visible, but the cosmic-ray tracks dominate the information content of the image. The scale is in 
ADU, and the image gain is 0.5 photoclcctrons per ADU. The image requires 3.1 bits per pixel to 
transmit after application of square-root compression and lossless bzip2 compression. 
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We compress our test images with a linear algorithm with code step A chosen to yield b = 1 



steps per noise a at the sky level of the images. As noted e.g. in Pence et al. (2009), most pixels 
are at sky level so a fixed b at sky level should yield similar compression performance. Indeed for 
images with no cosmic rays, we find that linear and square-root compression yield indistinguishable 
bpp values after the lossless step. As expected, the square-root algorithm performs better if the 
number of CR-impacted, high-signal pixels is substantial. In our nominal 200s exposure, the 
square-root algorithm uses 0.12 fewer bpp than decimation. For a 600s CR dose, the square-root 
algorithm requires 0.38 fewer bpp than decimation. This suggests that the square-root algorithm 
removes about 2 bits of noise from the average CR-impacted pixel. If an algorithm to remove 
cosmic rays before compression is implemented, then there would be little difference in compression 
ratio between the square-root and linear codecs. However such an algorithm would require more 
processing power than the lookup table of the square-root codec, and one would have to insure that 
objects of interest are not discarded with the cosmic rays. 

We should expect the square-root algorithm to gain advantage also in the case of fields with 
many objects brighter than the sky background, e.g. fields at low Galactic latitude or with large 
nearby galaxies. 

The square-root algorithm has practical advantages as well, in that it will retain the desired 
6 = 1 value as the sky background varies. The linear codec would need to change its step size A to 
accommodate different sky levels, e.g. from changing filters or zodiacal background. An on-board 
algorithm could determine the sky noise level and choose A appropriately for the linear codec. The 
square-root codec requires adjustment only when the bias level or gain of the detector ADC output 
change. 



7. Conclusion 

Future dark-energy surveys such as JDEM or E'ltc/icQ propose to survey a substantial fraction 
of the full sky with space-based visible/NIR imaging near the diffraction limit. A data compression 
scheme is needed to satisfy telemetry constraints, but the weak gravitational lensing measurements 
place very strict requirements on the degree to which the compression algorithm can bias the 
measured shapes or fluxes of galaxies. For this reason we choose to investigate lossy compression 
processes that operate on a pixel- by-pixel basis, fearing than algorithms examining pixel groups 
might induce correlated codec errors between pixels which bias the galaxy ellipticities. 

Square-root pre-processing is known to be a good means for compressing shot-noise-dominated 
data streams while keeping a fixed number h of bits per a of noise across the full detector dynamic 
range. The reduction of entropy by reducing the number of noise bits is essential to efficient subse- 
quent lossless compression algorithms. We derive an expression for bias in the naive reconstruction 



' http : //sci .esa. int/euclid 
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of square-root-compressed data and present a correction scheme, which we numericahy verify is 
successful to high accuracy. 

Since the noise induced by roundoff in the compression process increases the required observing 
time by 1/125^, we confine our examination to 6 > 1. We verify that fiux and shape measurements 
of simple galaxy models are biased by < 10""^ for a range of galaxy sizes, ellipticities, orientations, 
and signal-to-noise ratios. This is well within the weak lensing requirements. We verify this with 



simple postage-stamp galaxy images; Vanderveld et al. (2009) bound the shape biases determined 



from realistic simulations of the true extragalactic sky. 

At 6 = 1 the Shannon entropy of the square-root-compressed images varies from 2.1 bits per 
pixel for pure noise to 3.6 bits per pixel for images with the expected astronomical scene and 
triple the expected number of cosmic rays. Lossless compression algorithms such as bzip2 realize 
compress levels within 0.3-0.4 bits of the Shannon limit. For simulated images best resembling 
expected JDEM CCD data, we find compression to 3.1 bits per pixel using bzip2. The size of 
the transmitted data stream depends primarily on the degree of cosmic-ray contamination once 
6 = 1 is set: schemes to more efficiently compress or mask cosmic rays before compression might 
lead to lower bit rates. Higher cosmic-ray rates or more crowded astronomical fields, e.g. near the 
Galactic plane, will require more bits. Since pure noise images require 2.4 bits per pixel with bzip2 
at 6 = 1, we can state that compression to < 2.4 bits per pixel will require 5 < 1, and hence higher 
observing-time penalty. 

We thus validate square-root lossy compression as useful for JDEM data and capable of meet- 
ing requirements for low bias. Other approaches would be feasible as well. For instance simple 
decimation to a fixed code step A that equals the sky noise level should produce similar noise and 
bias properties. The presence of cosmic rays and/or bright objects in > 10% of the pixels gives 
the square-root codec a useful advantage in mean number of required telemetry bits per pixel. The 
square-root algorithm also has a robust and exceptionally simple implementation — a single lookup 
table — meaning that there is no need for on-board processing to recognize cosmic rays or determine 
the sky noise level. 

It is thus expected that visible/NIR dark-energy extragalactic imaging survey data can be 
transmitted to the ground with 2.5-4 bits per pixel, depending upon the cosmic-ray fiux level 
and the sophistication of onboard cosmic-ray identification or removal. This compression can be 
achieved with no significant bias of observed quantities, and with 8% penalty in observing time to 
achieve fixed signal-to-noise ratio. 
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